Worm algorithms for classical statistical models 
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We show that high-temperature expansions may serve as a basis for the novel approach to 
efficient Monte Carlo simulations. "Worm" algorithms utilize the idea of updating closed path 
configurations (produced by high-temperature expansions) through the motion of end points of a 
disconnected path. An amazing result is that local, Metropolis-type schemes may have dynamical 
critical exponents close to zero (i.e., their efficiency is comparable to the best cluster methods). We 
demonstrate this by calculating finite size scaling of the autocorrelation time for various universality 
classes. 
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PACS numbers; 75.40.Mg, 75.10.Hk, 64.60.Ht 

Metropolis scheme Q is usually the most universal and 
easy to program approach to Monte Carlo simulations. 
However, its advantages are virtually canceled out near 
phase transition points. It is believed that any scheme 
based on local g Metropolis-type updates connecting 
system configurations into Markovian chain is inefficient 
at the transition point because its autocorrelation time, 
T, scales as L^, where L is the system linear dimension 
and z is the dynamical critical exponent which is close to 
2 in most systems 

An enormous acceleration of simulations at the criti- 
cal point has been achieved with the invention of clus- 
ter algorithms by Swendsen and Wang However, the 
original method and its developments (both classical and 
quantum) ||^-^ are essentially non-local schemes, and we 
are not aware of any exception from this rule. 

In this Letter we propose a method which essentially 
eliminates the critical slowing down problem and yet re- 
mains local. The corner stone of our approach is the 
possibility to introduce the configuration space of closed 
paths. Closed-path (CP) configurations may be then 
sampled very efficiently using Worm algorithm (WA) in- 
troduced in Ref. for quantum statistical models in 
which closed trajectories naturally arise from imaginary- 
time evolution of world lines. In classical models the CP 
representation derives from high-temperature expansions 
for a broad class of lattice models (see, e.g. Ref. [0). 
In 2D, another family of WA may be introduced by con- 
sidering domain- wall boundaries as paths. 

We note, that our approach is based on principles 
which differ radically from cluster methods and, most 
probably, has another range of applicability. For one 
thing, the CP representation is most suitable for the 
study of superfiuid models by having direct Monte Carlo 
estimators for the superfiuid stiffness (through the his- 
togram of winding numbers [ pl] | ) which arc not available 
in the standard site representation. 

In what follows we first recall how high-temperature 
expansions work by employing Ising model as an exam- 
ple (still, trying to keep notations as general as possible). 
We then explain how WA is used to update the path 
configuration space. Next, we discuss specific implemen- 



tations of WA for XY-, and q = 3 Potts models, 

and comment on the special property of 2D models which 
allows an alternative CP parameterization of the configu- 
ration space. The efficiency of the new method is studied 
by looking at autocorrelation properties for six different 
universality classes. It is found that for 2D and 3D Ising 
models, 2D and 3D XF-models, and Gaussian model the 
t{L) scaling is consistent with the law t{L) — To+cln{L), 
i.e., its critical exponent is close to zero. For the two- 
dimensional q — 3 Potts model our data are consistent 
with the power law with z « 0.55 which means that the 
Li-Sokal bound z > a/v {a and v are the specific heat 
and correlation radius critical exponents) derived for the 
Swendsen- Wang algorithm is seemingly applicable to our 
method as well. 

Since high-temperature expansions for various models 
can be found in standard texts (see, e.g., Ref. |l^) we 
briefiy remind the procedure for the Ising model 
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where (3 — J/T is the dimensionless nearest- neighbor 
coupling parameter between spin variables Sj = ±1 and 
index b ~< ij > refers to the simple cubic/square lattice 
bonds (we will also use another notation: b = (i,i^), in 
which v enumerates bonds containing site i). Since H 
is additive, the corresponding Gibbs exponent factories 
in terms of exponents for each bond. Expanding each 
exponent in Taylor series allows one to perform summa- 
tion over site variables and to arrive at an expansion in 
powers of (3. The partition function, for example, takes 
the form 
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Note that summation includes all possible CPs, both con- 
nected and disconnected, with self-intersections and over- 
laps. Graphically, each elementary factor in the order- 
by-order expansion over the bond Hamiltonian can be 
represented by a line ascribed to the corresponding bond 
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|lO| , and the closed path constraint is imposed by sym- 
metry: the sum s^* is non-zero only if the number 
of bond lines connected to the site i is even. For future 
purposes it is convenient to denote the number of bond 
lines A'^!, = 0, 1, 2, . . . and call it the "bond state" . Then 
the "site state" is ki — J^u^^^i^^ ^^'^ configuration 
weight, Wcp, is given by 



Wcp 




(3) 



[Q{k) = 2 for the Ising model.] 

The spin-spin correlation function G{ii — 12) = 
{si^Si^) (we understand site indices as vectors) by def- 
inition equals to g{ii — 12) /Z where g{ii — 12) — 
^{a } ^ii^i2 ■ Configurations for g differ from 

those contributing to the partition function in only one 
respect: there are two special sites ii and 12 with odd 
number of bond lines connected to them, see Fig. |l|. The 
corresponding site states are given by fc^ = 1 -I- Ni^y. 
These points are the only places where the path may be 
disconnected. We will abbreviate the configuration space 
for g{i) as CPg, and note, that weights Wcp^ are given 
by the same expression (0). 



FIG. 1. A typical configuration with the "Worm". The two 
circles correspond to the extra Si^ and Si^ variables, and the 
solid line width is proportional to the bond state number A^;, 
(number of elementary lines); 0-th order terms are shown by 
dashed lines. 

The above consideration is actually a perfect setup for 
the Monte Carlo simulation since configurations for Z 
and for 5(0) have identical bond elements (for the Ising 
model they even have equal weights, i.e. g{0) = Z). 
Thus, if an ergodic Metropolis process is sampling g- 
configurations according to their weights, then the Monte 
Carlo estimators for g{i) and Z are just 
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The spin susceptibility is given by x = '^i9{^)/Z, and 
energy may be computed in two ways: (i) as the nearest- 
neighbor sum E = {L'^/2)(3J2\i\=i9{''')/Z, and (ii) using 
direct estimator, E = £/Z, which is non-zero only for 
the CP configurations contributing to Z: 



5(CP) =/?^iVb 
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Our results for r were obtained using this estimator. 

The WA emerges as an idea to update ^-configurations 
through the space motion of the end points ii and 12 only 
Q. The algorithm itself is extremely simple and con- 
sists of two elementary updates which are proposed in 
the context of a given configuration: if ii =12, then with 
probability po we suggest to "move" both end points to 
site io selected at random among all L"^ lattice cites, and 
with probability pi — I ^ po to "shift" the end point ii 
to one of the neighboring sites by selecting the direction 
of the shift at random; if ii ^ 12 we always suggest to 
"shift" ii. 

For the Ising model the move-update is automatically 
self-balanced and its acceptance ratio is unity, otherwise 
it would involve the ratio of Q-functions 



Prr 



Q{k, + 2) Qjk, - 2) 



Qih 



(6) 



In contrast, shift-updates also increase/decrease the bond 
state by "drawing/erasing" bond lines, and we select 
which way to proceed with probability 1/2. The shifts 
i ^ j with Nb^^ijy ^ Nb + I are balanced by shifts 
j ^ i with Nb ^ Nb — 1. By standard rules |l| we have 
acceptance ratios as 



Psi,ii^ j,Nb-^ Nb + 1) 
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Psh(*-^j,Wfc^7Vf,-l)=r-^ 



{Nb + l) Q{kj) 
Nb Q{h - 2) 



where 



p Qih) 



l/(2pi) if = ^2 = « , 

2pi if ii = i, i2 = j 
1 otherwise 
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The Ising model has a special property which can be 
used in practical applications to enhance the efficiency by 
truncating the configuration space. Namely, using iden- 
tity e'^^'^J = cosh/3[l + ta.nh PsiSj] one can restrict the 
bond summation to include terms Nb — 0,1 only. The 
corresponding modifications of the scheme are straight- 
forward. 

We now turn to the j'i/'l^-model described by the Hamil- 
tonian (ipi is a complex variable) 
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the two limiting cases of which are the Gaussian model, 
U — 0, and the Xy-model, U = 00. The procedure 
of factorizing the partition function expansion in terms 
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of the bond Hamiltonian is exactly as before. The only 
new ingredient which is not present in the Ising model, is 
that for each bond we have two different terms: ipiipj and 
ip*tpj, and therefore the expansion has to be performed 
for each of them separately. Graphically, this can be cap- 
tured by drawing lines with arrows, thus specifying each 
term by the arrow direction. Correspondingly, the bond 
state is defined by two numbers {Nj^^\ Nj^"^^) which tell 
how many lines and in which direction go along this bond. 

Since integrals / d^/'e'^l'^l'-'^l'^l = <5™^™'Q(fc = 

2m) are non-zero only for m — m' , we conclude that the 
configuration space for Z is defined on CP and site states 



are given by h = (^i'J + ^i'J) =even. [The Q{k) 
integrals are easily tabulated prior to the simulation.] 

In close analogy with the previous case, the con- 
figuration space for the two-point correlation function 
9{ii -12)= J llidipi i^iii^l^ e"-^/^, is given by CPg con- 
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taining two special sites with ki — l + Ylv {j^t 
WA updates consist of "moves" which change the loca- 
tion of end points when 11—12 and "shifts" of ii which 
increase/decrease numbers n[^\ and N^^\. Acceptance 
ratios are given by Eqs. (^[^) where iVb is understood as 
one of the numbers describing the bond state. An estima- 
tor for the winding number (which has exactly the same 
meaning as in worldlinc quantum Monte Carlo jll]]) of 
the CP configuration is given by the difference between 
the bond numbers 
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where a = x, y, z, . . ., and notation ba is used to specify 
bonds connecting sites in direction a. 

The q — 2> Potts model (which we considered to com- 
pare with the exhaustive cluster method study by Li and 
Sokal is described by 



note that in lattice models with discrete site variables one 
can unambiguously specify the state of the system (up to 
symmetry transformations) by drawing domain bound- 
aries, which, in 2D, may be considered as a CP config- 
uration. For the Ising model the configuration weight is 
simply Wcp = J^^^ e'^^^^''"^^ where the bond state takes 
values Nb = 0, 1. 

To implement WA for efficient sampling of domain 
boundaries we formally enlarge the configuration space 
to include CPg configurations with two end points and 
proceed exactly as discussed above. The only difference 
is that now "open" boundaries with ii ^ 12 have no phys- 
ical meaning and serve just for updating purposes. Our 
results for the 2D Ising model in Fig. g were obtained 
using domain-wall representation. 
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2D Ising model 



2 3 4 5 6 ln(L) 



2D Xy-model 
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Introducing phase variable kpi — (27r/3)si we can rewrite 
Eq. (O) (up to a constant term) as 



f 



P' (e*^'^'-'^^-) + c.c.) , P' = /5/3 . (13) 
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3D Gaussian model 
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q=3 Potts model 



3 4 5 ln(L) 



With this form one may immediately proceed with the 
WA along the fines described above for the |V'|'*-niodel. 
However, more efficient scheme results from the iden- 
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2/3' cos(c^i— c^j) 
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7COs((/3i — ipj)], where A = 
{e^P' + 2e-P')/Z and 7 = 2{e^'^' - e-P')/iA, which al- 
lows to restrict summation over bond states to just three 
values Nb = (0,0), (1,0), (0,1). 

We would like to comment that high-temperature ex- 
pansions are not the only procedure to arrive at the CP 
configuration space and WA. To underline this point we 



FIG. 2. Autocorrelation times for various universality 
classes. The 3D Ising model is fitted to r = —4.3 -I- 9.2 ln(L), 
and z{L = 64) = 0.18 (see text). The 2D Ising model is 
fitted to r = -7.2 + 6.21n(L), and z{L = 512) = 0.25. 
The 3D J^y-model is fitted to r = 1.7 + 3.85 ln(L), and 
z{L = 128) = 0.2. The 2D XY-model is fitted to 
T = 1.85 -f 2.05 ln(L), and z(L = 640) = 0.16. The 
3D Gaussian model is fitted to r = 18.9 + 5.81n(L), and 
z{L = 64) = 0.17. The g = 3 Potts model in 2D is fitted 
to the power law r = 4.3L'^'^^. 
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In Fig. g we present data for the energy autocorrelation 
time. In all cases, except g = 3 Potts model in 2D, they 
are consistent with the linear in In L behavior. Of course, 
we may not exclude that logarithmic scaling is actually 
an intermediate length-scale behavior and at some larger 
L it crosses over to the power law, t{L) ^ U- , with small 
z. To get an estimate (upper bound) for the possible 
dynamic exponent we mention in the figure caption the 
slope z{V) = dln(r)/ciln(L) at the largest simulated L. 
For the (7 = 3 Potts model our data scale as r ~ /^0-55 
for L > 64, although the dynamic exponent is showing 
a systematic decrease with L. This has to be compared 
with the value z = 0.515 for the Swendsen-Wong algo- 
rithm . It seems that the Li-Sokal lower bound for the 
dynamic exponent z > a/v also applies to our method. 
For all models discussed here a/v is very small with ex- 
ception of the f? = 3 Potts model were a/v = 0.4. 

To summarize, we found that WA Monte Carlo 
schemes working with closed path representations of clas- 
sical statistical models are as efficient as the best clus- 
ter methods in reducing the problem of critical slowing 
down near transition points. The new approach is gen- 
eral enough to find applications other than discussed in 
this Letter, and may be used to study quantities like su- 
perfluid density for which other methods may not have 
direct estimators. 
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